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The renormalization of electronic eigenenergies due to electron-phonon coupling (temperature dependence 
and zero-point motion effect) is sizeable in many materials with light atoms. This effect, often neglected 
in ab-initio calculations, can be computed using the perturbation-based Allen-Heine-Cardona theory in the 
adiabatic or non-adiabatic harmonic approximation. After a short description of the numerous recent pro¬ 
gresses in this field, and a brief overview of the theory, we focus on the issue of phonon wavevector sampling 
convergence, until now poorly understood. Indeed, the renormalization is obtained numerically through a 
(usually slowly converging) q-point sampling inside the Brillouin Zone. For q-points close to T, we show that 
a divergence due to non-zero Born effective charge appears in the electron-phonon matrix elements, leading 
to a divergence of the integral over the (phonon) Brillouin zone for band extrema. Although it should vanish 
for non-polar materials, unphysical residual Born effective charges are usually present in ab-initio calcula¬ 
tions. Here, we propose a solution that improves the coupled (electronic) k-point convergence dramatically. 
For polar materials, the problem is more severe: the divergence of the integral does not disappear in the 
adiabatic harmonic approximation, but only in the non-adiabatic harmonic approximation. In all cases, we 
study in detail the convergence behavior of the renormalization as the q-point sampling goes to infinity and 
the imaginary broadening broadening parameter goes to zero. This allows extrapolation, thus enabling a 
systematic way to converge the renormalization for both polar and non-polar materials. Finally, the adi¬ 
abatic and non-adiabatic theory, with corrections for the divergence problem, are applied to the study of 
five semiconductors and insulators: a-AlN, /I-AIN, BN, diamond and silicon. For these five materials, we 
present the zero-point renormalization, temperature dependence, phonon-induced lifetime broadening and the 
renormalized electronic bandstructure. 


I. INTRODUCTION 

The theoretical understanding of the effects of the 
electron-phonon coupling on the electronic structure and 
the capability to compute them have a long and chaotic 
history that started in the early fifties. Over the years, 
these effects have been computed using three types of 
methods, with different advantages and drawbacks: (1) 
as a time average of the bandgap using first-principles 
molecular dynamics (MD) simulations; (2) through the 
frozen-phonon (FP) method, which weights the eigenen- 
ergy change along the phonon modes with a Bose- 
Einstein distribution; (3) thanks to the diagrammatic 
method of many-body perturbation theory. For a histor¬ 
ical review, the reader can consult Ref. [TJ in which these 
three types of methods are compared to each other, at 
the harmonic level. 

In the present contributiom we rely on the Allen- 
Heine-Cardona (AHC) theorySSl to compute the zero- 
point motion renormalization as well as the temperature 
dependence of electronic eigenenergies. The AHC the- 
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ory originates from the diagrammatic method of many- 
body perturbation theory. It has been applied in sev¬ 
eral recent milestone contributions in the field, includ¬ 
ing the computation of temperature-dependence of the 
optical propertied, the computation of the surprisingly 
large zero-point renormalization (ZPR) of the diamond 
bandgapl^, the demonstration of large non-rigid ion cor¬ 
rections for moleculed, the inclusion of dynamical effects 
beyond the adiabatic approximatiorPtH^, the study of the 
anharmonic electron-phonon contribution to the indirect 
bandgap of diamoncP^, and the inclusion of electronic 
many-body effects (in the GW approximation) in dia- 
moncP^, noticing a large increase of the renormalization 
with respect to density-functional theory (DFT). Also, 
we think that the confusion in the theoretical under¬ 
standing of the relationship between MD, FP, and AHC 
as well as the inaccuracies in first-principles software im¬ 
plementations of AHC have been largely eliminated in 
two recent publication^^. 

One of the major issues when performing AHC calcu¬ 
lations is the slow convergence with respect to phonon 
wavevector sampling of the Brillouin Zone (BZ)P, refer¬ 
eed to as q-point sampling from now on. To accelerate 
this convergence, a small imaginary component iS (which 
can be inferred as a finite lifetime for the unoccupied 
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electronic states due to thermal effects) is often used. 
However, this imaginary parameter is ad hoc rather than 
ah-initio. Also, the convergence problem is even more 
severe with the MD and FP methods, as supercells have 
to be used to sample the phonons wavevectors, thus dra¬ 
matically increasing the computational time and mem¬ 
ory required. Actually, numerical convergence for the 
MD and FP methods cannot really be reached in three- 
dimensional solids, in contrast with finite systemj ^ff^ff^ l 

In this paper, we highlight that when trying to con¬ 
verge the ZPR with respect to q-point sampling for van¬ 
ishing iS in the AHC simulations, the ZPR diverges. For 
non-polar materials, such unphysical divergence is at¬ 
tributed to a residual Born effective charge, which stems 
from the finite k-point sampling. We propose a solu¬ 
tion to this problem and devise a systematic way to con¬ 
verge the ZPR for vanishing iS. For polar materials, the 
problem is more profound. Indeed, the divergence in the 
adiabatic AHC approach is not simply numerical, but 
indicates a breakdown of the AHC approach. A similar 
problem should also be present in the MD and FT meth¬ 
ods. On the other hand, the non-adiabatic AHC theory 
naturally leads to non-diverging quantities. 

This paper is organized as follow. First, a short re¬ 
minder of the AHC theory is presented in section |TTj In 
Sec. |HH the bottleneck of the q-point convergence is dis¬ 
cussed, the divergence problem of the ZPR at large q- 
point density is explored and a solution is proposed. We 
also device in sections |IV] and |V] a systematic and pa¬ 
rameter free way to extrapolate the ZPR (without i6). 
Finally in section [Vij we present the temperature depen¬ 
dences, the zero-point motion renormalizations, as well 
as the phonon-induced lifetimes for five semiconductors: 
Of-AIN, /3-AlN, BN, diamond and silicon. 


II. REVIEW OF THE ALLEN-HEINE-CARDONA 
FORMALISM 

A. The AHC theory within the adiabatic harmonic 
approximation 


The temperature-dependent renormalization of the 
electronic eigenenergy Snk. for band n and wavevector k 
can be written in the adiabatic harmonic approximation 
as a sum over the BZ of the phonon contributions for 
each wavevector 
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where m is the phonon branch, T is the temperature, Nq 
is the number of wavevectors used to sample the BZ, 
is the phonon frequency, nmq{T) = — is the Bose- 

e '‘bT _i 

Einstein distribution, Um,Kai^) is the eigendisplacement 
vector of atom k in direction a associated to the phonon 
mode, and d/dRi^a is the derivative of a quantity with 
respect to the displacement of atom k of the unit cell I 
in the direction a. 

The quantity Ae„k(F) A £nk{T) — Enkp] is the dif¬ 
ference between the temperature-dependent eigenenergy 
and the eigenenergy at ground-state atomic positions. 
We distinguish the temperature-dependent e„k from the 
atomic-position-dependent Snk by using () in the former 
and [] in the latter. The difference between £nk{T = 0) 
and enk[0] is called the zero-point motion renormaliza¬ 
tion (ZPR). 

In a mean field approximation like the Density Func¬ 
tional Theory (DFT), the eigenenergies are the expecta¬ 
tion values of the Hamiltonian TFk.k of the system 
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with the periodic part of the electronic wavefunc- 
tions. Using perturbation theory to obtain the second- 
order derivative with respect to atomic displacements of 
such eigenenergies, Eq. ([^ can be rewritten as 
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where (c.c.) stands for the complex conjugate of the pre¬ 
vious terms within parenthesis (), and where we use 
the following notation for the derivative with respect to 
atomic positions of an arbitrary quantity X 
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NbvK being the number of primitive cells of the peri¬ 
odic system defined by the Born-von Karman boundary 
condition^i^. 

The first term within {} in Eq. Q is called the Debye- 
Waller (DW) term 
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while the remainder constitutes the Fan term 


with 


TK.a (q) = - 

Ki 7 Z 


du. 


dHu] 


iidR^aiq) 9^K'7(q) 

+{Ka) O (k' 7)) + (c-c.) 


“nk / 


(7) 


The change of eigenenergy due to a specific phonon 
mode (e.g. Eq. @) thus becomes 
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with the Fan contribution given by 
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and the Debye-Waller contribution given by 
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At this point, no approximations beyond the adiabatic 
and harmonic ones were made. However, the calculation 
of the Debye-Waller term (Eq. Q) requires the second- 
order derivative of the Hamiltonian, which is a computa¬ 
tional bottleneck within the density functional perturba¬ 
tion theory (DEPT) framework. To overcome this issue, 
we make the rigid-ion approximation (RIA) as is usual 
within the AHC theory B®!. We begin by splitting the 
DW term into two parts. The first part contains all con¬ 
tributions that can be computed using first-order deriva¬ 
tives of the Hamiltonians, while the second part contains 
the remaining contributions 
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Moreover, in our calculations, all Ean-like contribu¬ 
tions are obtained within DEPT and can thus be written 
as follow 
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where the summation over energetic bands (above M) 
has been replaced by the solution | Ak+q ) of a 

linear equation as proposed by Sternheimeili^ and applied 
to this problem in Ref.[Sl The definition for the projector 
Pck+q and active space M as well as the description of the 
linear equation to be solved can be found in the appendix 
of Ref. [fl 

We finally obtain the adiabatic temperature-dependent 
renormalization in the RIA by neglecting the non-RIA 
contribution as defined by Eq. (13), which yields 
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where a small imaginary component i5 is usually intro¬ 
duced in the AHC equation to smooth the energy de¬ 
nominators. For example, in the case of diamond, several 
authors have used an i5 of 100 meV to accou nt for the 
finite lifetimes of the electronic statej^ ^^ l ^^ l However, 
the theory must also be valid (apart from controlled nu¬ 
merical instabilities) for vanishing i5. This point will be 
further discussed in section Ell 

B. Beyond the Rayleigh-Schrodinger perturbation theory 

Phonons alter the one-electron energy bands e„k in 
two ways: there is a shift Asnk and a lifetime broad¬ 
ening llTnk- As seen in the previous sub-section, the 
adiabatic approximation leads to a real renormalization 
of the eigenstates. The study of the lifetime broadening 
requires an extension of the adiabatic theory. 

In 1978, Allen generalized his earlier worlP derived 


within the standard Rayleigh-Schrodinger perturbation 
theory to include finite phonon frequencies using many- 
body perturbation techniqueJi^. These techniques de¬ 
scribe excitations in terms of spectral functionJiSl^ where 
quasiparticules cannot always be unambiguously identi¬ 
fied, with the associated well defined eigenenergies. In 
this work, following AllerP^, rather than obtaining the 
full spectral function to describe the electronic excitation, 
we suppose that their description in terms of quasipar¬ 
ticles is still valid and evaluate the associated eigenen¬ 
ergies by correcting the DFT eigenvalues to first-order 
in perturbation theory, taking the self-energy evaluated 
at w = as the perturbation. Complex eigenener¬ 
gies are obtained within this generalization, that we refer 
to as the “non-adiabatic” extension of the AHC theory 
£„k(r,a;:pll2l. 

We therefore obtain the following equation, based on 
electron-phonon matrix elements already calculated for 
the adiabatic renormalization 
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where /x is the chemical potential, /„k is the electronic 
occupation of the wavevector k at band n and where a 
convergence study on M is required for the Fan term 
due to the fact that the Sternheimer solution neglects 
the phonon frequency ujmq while the sum over the active 
space does not. 


imaginary part of the complex Fan self-energy 
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The phonon-induced lifetime broadening I/r^k is the 
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where <5 is the Dirac delta (broadened for numerical rea¬ 
sons). 

The phonon-induced lifetime broadening in the adia- 
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The adiabatic and non-adiabatic renormalizations, 
Eqs. IT^ and (16), as well as the adiabatic lifetime 
Eq. (18) have been coded in the ABINIT software 
(vT-lljrand will be used in the following sections. 


III. PHONON WAVEVECTOR SAMPLING AND THE 
DIVERGENCE PROBLEM 


However, in non-polar materials, there are theoretically 
no such effective charges and there should therefore be 
no divergence of type (ii). Since divergences of type (i) 
have a finite integral when no divergence of type (ii) are 
present, the q-point sum should convergence to a finite 
value for non-polar materials. 

In practice, however, we observed a non-physical di¬ 
vergence of the ZPR for large q-points densities. This 
effect can be clearly seen in Figure where the q-point 
density dependence of the adiabatic T^g ZPR of diamond 
(calculated using Eq. ( [T^ ) for vanishing iS exhibits a di¬ 
vergent behavior when (5 is 1 meV or 0.01 meV, hardly 
seen for 6 equal to 50 meV or 100 meV. 



1/^, 


A. Potential breakdown of perturbation theory 


Quantum mechanical perturbation theory can break¬ 
down when vanishing denominators appear in the per¬ 
turbation series. This can happen in the present case, as 
the short-hand form of Eq. (15) is 


A (adiabatic,RIA)/rTAX \ ' l^-^-^(*l)l 

A^„k (^)°^E7 o)_jor 


(19) 


Actually, there are two types of potential divergences 
in Eq. (19): (i) when and (ii) when 

the electron-phonon matrix elements GKK{q) diverge, 
which happens when the sum of Born effective charges 
does not vanish, as we shall see {GKK{q) is then pro¬ 
portional to ^). 

In practical calculations, the q = 0 contribution from 
the same band (the denominator being thus zero) is not 
included in the summation. Also, in case of degeneracies, 
the terms with zero denominators are ignored. However, 
the integral of these divergences still needs to be obtained 
through the q-point summation. For this reason, the 
numerical convergence of the adiabatic ZPR of diamond 
with respect to q-point density is slow and requires large 
q-point grid^. This problem is often assessed in practice 
by adding an ad-hoc i6 to the denominator of Eq. (19). 

The dipoles present in polar materials induce a Born 
effective charge, which describes the coupling between 
the electric field generated by the dipoles and the ionic 
motion. Such Born effective charges lead to a ^ behavior 


of the electron-phonon matrix elements (GKK). Diver¬ 
gences of type (ii) are therefore present in these materials. 


FIG. 1. (color online) Adiabatic ZPR of the state of dia¬ 
mond with respect to the q-point grid density for decreasing 
values of iS before restoration of the charge neutrality. The 
size of the q-point grid is Nq x Nq x Nq. 


Actually, this divergence is attributed to a residual 
electric field connected to the breaking of the Born ef¬ 
fective charge neutrality sum rule in non-polar periodic 
solids. This leads to the presence of divergences of type 
(ii), which combined to these of type (i), give an infi¬ 
nite integral around q = 0 , thus making the q-point sum 
diverge. 

This residual electric field is due to the finite k-point 
grid used within DFT. Indeed, the first-order density is 
obtained thanks to a discretized integral on the BZ (see 
Eq. (B9) of Ref. EOl) . This first-order density in turn de¬ 
termines the electric field and the Born effective charges 
(see Eq. (42) of Ref. HD). Such residual electric field that 
breaks the charge neutrality is found to converge to zero 
exponentially but is nonetheless substantial at k-point 
grid usually sufficient to converge other relevant quan¬ 
tities. An example of this slow convergence is given in 
Table |T] for the case of diamond (a non-polar material). 


B. Restoration of the charge neutrality 

In this section, we present a scheme to numerically re¬ 
move this spurious electric field, and thus considerably 
speed up the ZPR convergence with respect to the elec¬ 
tronic wavevector sampling. To this end, let us study the 
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Number of k-points 

Born effective charge 

8 

-2.5406 

64 

-0.3514 

216 

-0.0534 

512 

-0.0080 

1000 

-0.0011 


TABLE I. Born effective charge of one carbon atom in dia¬ 
mond, for different electronic wavevector samplings. 


impact of a small Born effective charge on the matrix el- 
ements of present in Eq. ( |15[ ). We introduce 

the following short-hand notation for derivatives, which 
matches the one used in Ref. [20] 


_A ttW 
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As mentioned in Ref. [201 within the pseudopotential 
framework, can be decomposed into a first-order 

change of the non-local, local, Hartree and exchange- 
correlation potentials: 
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where the bar symbol above a quantity X means that it 
is the periodic part of X 

X(i)(r)=e-»'i-X(i)(r), (22) 

and where both ^ and ^ diverge as with oppo¬ 
site signs. 

To make this more explicit, we express ^ as 


^/oc q(G) = ^(G + q)ae-*(°+'>) -^^r(G -b q), (23) 

where fig is the volume of the unperturbed unit cell, 
the vector position of the atom k in the unit cell and with 


^^(q^0) =-^+C. + 0(q^). (24) 
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To also explicit the same behavior in we express it 
as 


.y(l) 


dl) 
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(25) 


is presented in Appendix [S] of this paper and can be seen 
in Eq. ( A77| ) in the G = 0 limit. 

Using this knowledge, we can renormalize the Hartree 
term as follow (see Eq. (A81)) 


-ren(l) 

(0) 
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where Z* is the Born effective charge, is the macro¬ 
scopic static dielectric tensor for the electronic system 
(where the ions are considered fixed) and Za^g is the av¬ 
eraged Born effective charge 


-'a/3 
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N^t 
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(27) 


where Nat is the number of atoms in the primitive cell. If 
the charge neutrality sum rule was fulfilled, the averaged 
Born effective charge should be exactly zero. 

With this renormalization, the term correctly 

cancels the v^l ^ when q —0. Figure 2 clearly shows the 
faster convergence rate of the ZPR with respect to the 
density of the k-point grid obtained with this renormal¬ 
ization for the specific case of the first band of diamond 
at k = L. To highlight the divergent behavior of the 
ZPR with respect to q-point grid density without the as¬ 
sociated high computational cost for q-point integration, 
only 6 symmetry equivalent q-points are used in the sum. 
The q-points are chosen close enough to zero to show the 
divergence: q = 



where riq^ cx |q| when q —> 0 . 

The vlll q(G) of Eq. (23) has an explicit algebraic form 
for q —0, where is the number of valence electrons of 
the atom k described in the pseudopotential. Therefore, 
a residual electric charge can only affect the first-order 
density hq^ in Eq. (251. The derivation of the impact of 
a residual Born effective charge on the first-order density 


FIG. 2. (color online) Example of the contribution to the 
ZPR of 6 q-points with respect to the densification of the 
k-point grid of the first band of diamond at k = L. Only 
6 symmetry equivalent q-points were used to calculate 

the ZPR. The calculations were done with and without the 
renormalization of Eq. (|26|). 
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IV. BEHAVIOR OF THE q-POINT CONVERGENCE 

After enforcing the charge neutrality by application of 
Eq. ( [ 2 ^ , the theoretical rate of convergence of the ZPR 
can be analyzed when the number of q-points along a 
side of the Brillouin Zone Nq increases (the total number 
of q-points in the Brillouin Zone is N^)- After isolating 
the divergent behavior, we will analyze analytically, and 
numerically on simple models, the rate of convergence of 
the ZPR. 

We will observe that the q-point convergence can be ei¬ 
ther constant, linear {1/Nq) or divergent {Nq) depending 
on the state that is renormalized, the use of the adia¬ 
batic (Eq. ( |l^ ) or non-adiabatic (Eq. ( [T^ ) equation as 
well as the polar or non-polar nature of the material. The 
table (|n| gives a summary of those behavior. 


Non-polar 

Polar 


Cases 


VBM/CBM 

other 

VBM/CBM 

other 


q-convergence 

Adiabatic Non-adiabatic 


1/Nq lIVBir 

flat (|rvx||l 


NaWW 


1/Nq Ifrv^2l 


flat 

flat 


iIVAli 


lIVTO 


i/W | 1V 5^ 

l/Nq iTVM 


TABLE II. Convergence behavior with the densification of the 
q-grid and behavior for vanishing iS at converged or extrapo¬ 
lated q-grid. The only case that diverges is a polar material at 
the VMB/CBM nsing the adiabatic equation. The referenced 
sub-sections are given in parenthesis. 


A. Rapid convergence with N, 


The corrected (using Eq. (26 


l) electron-phonon matrix 
elements of non-polar materia s have no strong q-point 
dependence. Also, if the state of interest (nk) is a valence 
band maximum (VBM) or a conduction band minimum 
(CBM), the band dispersion is quadratic in reciprocal 
space around k, and therefore e„k — ^n'k+q behaves as 


1. Non-polar materials in the non-adiabatic approximation 
at VBM/CBM 


For a non-polar material within the non-adiabatic ap¬ 
proximation (Eq. im), we can model the q-point behav¬ 
ior of the ZPR of the VBM with 


lim 

(5-)-0 





^nk ^n'k+q T ^q T 


(28) 


parabolic behavior of the extrema leads to 



A Taylor expansion around Qc = 0 reveals that 


47r- 




(33) 


which means that the contribution from the integration 
around q = 0 is simply proportional to the volume of 
integration, as expected given the non-divergent nature 
of the integrand. 

Thus, neglecting the q = 0 contribution in the q-point 
sum of the non-adiabatic ZPR for a band extrema of 
a non-polar material causes an error proportional to 
whereas discretization of the q-point integration over the 
Brillouin Zone with the rectangle method causes an error 
proportional to Therefore, the error caused by the 


neglected q = 0 contribution is not visible in the global 
convergence behavior for this ZPR. This behavior will be 
referred to as “flat” convergence with respect to q-point 
grid density from now on. 

Additionally, we can numerically integrate Eq. (28) on 
a three dimensional grid of q-points using 


Nj) ^ q'^ +UJ + iS' 

1 q 


(34) 


where V is the volume of integration and where the ele¬ 
ment of volume is inversely proportional to the number 
of q-points Nq needed to discretized the grid. An exam¬ 


ple of the convergence of Eq. (34) with is shown on 
the top of Figure!^ for qc = 0.5 and uj = 0.01. 

This function converges very quickly with increasing q- 
sampling as can be seen e.g. on the center of Figure[6a|for 
the VBM of diamond using the non-adiabatic equation. 


2. Non-polar materials in the adiabatic or non-adiabatic 
approximation for a non-extremal band 


where q is integrated in a sphere of radius qc- The same 
derivation applies for the CBM with — Wq, the energy 
difference e„k — Cnk-i-q being negative. The phonon fre¬ 
quency shifts the poles of the function, so that the in¬ 
tegrand is analytic over the domain of integration. The 


If the state that we would like to renormalize is not 
a VBM nor a CBM, the denominator of the adiabatic 
(Eq. @) or non-adiabatic (Eq. (dU)) equations will 
be small when the state that we consider (e„k) has al¬ 
most the same energy as another state (en'k-i-q), minus 






































4.5 






iS 

FIG. 3. (color online) Behavior of the numerical integral of 
Eq. ( |34[ ) with respect to the q-point density and S for qc = 0.5 
and ui — 0.01. 


a phonon frequency (w^q) in the non-adiabatic frame¬ 
work. As a result, the integrand in Eq. (281 is not ana¬ 
lytic anymore in these cases and a non-zero imaginary iS 
is required to avoid numerical instabilities. 

To get a deeper understanding of the behavior of non¬ 
extremal points, we will model the energy difference 
(^nk — £n'k+q) by a shifted parabola with its minimum 
at qo 


Sk - ffk+q = (q- qo)^ - Qo- 


(35) 


The ZPR has poles on a sphere of radius qo centered on 
q = qo (chosen to be on the z-axis) and passing through 
the origin. The integral on the spherical shell between 
radii of values qo — A. and go + ^ gives a contribution 
that is linear with A (see Appendix B 1), as would any 
regular function when integrated over a spherical shell. 
This indicates that the integration of these poles will not 
contribute an error of higher order than the reminder of 
the numerical integration. Neglecting the q = 0 con¬ 
tribution in the numerical integration leads to an error 
proportional to ^ in the non-polar case (see Appendix 


B 2). This leads to a q-point grid convergence that is flat 


for the non-polar case in the adiabatic framework. 

In the non-adiabatic case, the w to be added to the 


right-hand side of Eq. (35) only slightly reduces the ra¬ 


dius of the sphere at qo (which thus does not touch 
the origin anymore) and the conclusion for the adiabatic 
case remains valid for the numerical integration over this 
sphere. The integrand at q = 0 becomes analytical in 
this case, so that the convergence behavior remains ef¬ 
fectively flat. Therefore, in practice, the q-point inte¬ 
gration required to evaluate the ZPR can be considered 
converged when the ZPR does not change significantly 
with denser grids. An example of this type of conver¬ 


gence is given at the top of Figure [5a| for diamond in the 
adiabatic framework at a non-extremal energy. 

In other cases, the discretized integral does not con¬ 
verge as quickly as the rectangle method for an analytical 
integrand. It sometimes converges linearly {oc (see 
subsection IV B) or even diverges (oc Nq) (see subsection 


IV C). 


B. Convergence proportional to the inverse of N, 

1. Non-polar materials in the adiabatic approximation at 
the VBM/CBM 

Non-polar materials with a parabolic energy dispersion 
(VBM or CBM) have a q-dependence for the adiabatic 
ZPR that behaves as 


Sq 


1 


q^ -\- iS 


= 4.9= « i 


(36) 


when iS = 0. Therefore, neglecting the q = 0 contribu¬ 
tion in the numerical integration yields an error propor¬ 
tional to that dominates the error of the rectangle 
method. We will call this type of convergence “linear” 
here. 

The rate of convergence with q-densification can be 
numerically tested by summing this function on a three 
dimensional grid of q-points 


— F 

Af3 






(37) 


where the q = 0 term has been omitted in the sum for 
numerical reasons (as the expression must stand for van¬ 
ishing S). The numerical integral of Eq. (371 is shown 
on Figure and converges towards 27r for qo — 0.5, as 
expected. 


2. Polar materials in the adiabatic (non-extremal point) 
and non-adiabatic approximation 


In the case of polar materials, the GKK behaves as 
1/q^ for small q, but this time, the divergence is physical 
and comes from the physical finite Born effective charges. 
In this case, at the VBM or CBM, the non-adiabatic 


Eq. (16) shifts the pole inside the bandgap. Therefore, 


the parabolic eigenenergy dispersion is not problematic 
anymore. This can be represented by 


5R 


9c 1 

j3 3 

0 q^(uJ + tS) 


(38) 


This case behaves as Eq. (37) and we have the same linear 


behavior because the diverging q = 0 term is removed 
from the numerical sum. 
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iS 


FIG. 4. (color online) The numerical integral of Eq. (371 
converges toward 2n as expected. We can see that the con¬ 
vergence in q points and iS is rather slow. 


In the case of points that are not at a VBM or a CBM, 
the integration of the ZPR in the adiabatic and non- 
adiabatic approximations exhibits a behavior that is lin¬ 
ear with 1/Nq. This is explained by the removed q = 0 
contribution and the integral is linear with the radius of 
the sphere (as shown in Appendix |B 2[ ), which is propor¬ 
tional to l/Nq. 


C. Increasing Nq does not lead to convergence 


In the case of polar materials, when we consider the 
renormalization of a state at the VBM or the CBM, the 
adiabatic equation diverges. Indeed, 


V y. 1 


(39) 


diverges as Nq. An example of this case is given in the 
middle left Figure [T6b| for the VBM of boron nitride using 
the adiabatic equation. 

More specifically, neglecting the q = 0 contribution in 
the numerical integration of Eq. (391 gives a convergence 
behavior that can be modeled by the spherical integration 
of the summand of Eq. (39) in a shell from qc/Nq to qc- 
This gives, for d = 0 



As l/Nq goes to 0, the value of this integral indeed di¬ 
verges linearly with the number of division Nq. This 
shows that, for polar materials, only the non-adiabatic 
equation can be safely used. 


V. BEHAVIOR OF THE iS CONVERGENCES 


After enforcing the charge neutrality by application of 
Eq. (26), the theoretical rate of convergence of the ZPR 


can be analyzed when the small imaginary parameter iS 
tends to 0. After isolating the divergent behavior, we will 
analyze analytically, and numerically on simple models, 
the rate of convergence of the ZPR. 

We will observe that the <5 convergence can be square- 
root, linear or Lorentzian-like, depending on the state 
that is renormalized, the use of the adiabatic or non- 
adiabatic equation as well as the fact that the material 


is polar or not. The table (III) gives a summary of those 
behavior. 


Non-polar 


Polar 


Cases 


VBM/CBM 
other 


q-convergence 


Non-adiabatic 


Lorentzmn Jy C 11 

Tf 


S ||VB^ 


VBM/CBM l/y JS ilVD t Lorent zian~W C2 1 
other <5 ifVB^ 5 ||VB^ 


TABLE III. Convergence behavior for vanishing i5 at con¬ 
verged or extrapolated q-grid. The only case that diverges is 
a polar material at the VMB/CBM using the adiabatic equa¬ 
tion. The referenced sub-sections are given in parenthesis. 


A. Convergence proportional to the square-root of 5 


For a non-polar material at the VBM or CBM, in the 
adiabatic framework, we can determine the ZPR depen¬ 
dence on the (finite) 5 value by analytically integrating 


Eq. (37) 


3 ? 


d^ 


q 


iS 


= drr 




6^ 



dq §2 
VS'^ + I 


(41) 


Carrying the integration and making a Taylor expansion 
for small y yields 

'i-Kqc + A'kC'/S + 0{5), (42) 

with C a constant. This result matches the behavior 
observed in our numerical integration in the bottom of 
Figure]^ 


B. Convergence proportional to 5 

1. Non-polar materials in the adiabatic and non-adiabatic 
approximation at a non-extremal point 

For non-polar materials, when considering a k-point 
other than the VBM or CBM, the adiabatic equation 
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can be modeled as 

3fJ 


d^q 


£(q) + i(5’ 


(43) 


where we have defined e(q) = Ek+q — £k- 

Since the small iS will only affect the integrand around 
the pole at e(q) = 0, we can determine the 6 dependence 
of this model ZPR by considering only a small range of 
energy rj around it. We can then re-write Eq. (43) as 


dEg{E)f{E), 


(44) 


with g{E) = ^ 5{E — e(q)) the density of states and 


f{E) = 


1 


E + i5' 


(45) 


We can make a Taylor expansion of the density of 
states around E = 0 (since rj is small) 

giE)=g{0)+g'i0)E + O{E^). (46) 


The even terms in Eq. (46) will not contribute because 
5ft/(E) is an odd function, therefore making no contribu¬ 
tion to the integral. Therefore, the leading term is the 
first-order one. The integral can thus be re-written as 


n dEg'{Q)^ 

L ^ S + 1 


25g'(0)( 


I—-‘(|))' 


(47) 


to 


Making a Taylor expansion of Eq. (47) for small ^ leads 


2. Polar materials in the adiabatic and non-adiabatic 
approximation at a non-extremal point 


For polar materials, we can use the shifted parabola 
model (defined in Eq. (35)) in Eq. (43), multiply the 
integrand by ^ and for the non-adiabatic framework, 
add w to i5. With in this model, it is possible to show 
(see Appendix |B 3| ) that, for small 6, the ZPR converges 
linearly with 5. In practice, the value that is reached 
by the linear regime is very high and tends to infinity at 
VBM or CBM. 


C. Convergence in 5 proportional to a Lorentzian 

1. Non-polar materials in the non-adiabatic approximation 
at tbe VBM/CBM 


When considering the non-adiabatic equation at the 
VBM or CBM of a non-polar material, the ZPR behavior 
with respect to 6 is 


5ft // / d\- 


— = 47r3fi / dq 


q^ +u: + id 

which gives 

47r5ft^gc — ^/oj + j(5tan“^ ^ 


Q q^ +0J + id 


— if dc 


y/uj + id 


(54) 


(55) 


2g'{0)1] - dg'{0)71, (48) 

where we can see that this function is linear in d. 

In the framework of the non-adiabatic theory, the in¬ 
tegral is changed to 

J dEf{E + Lu)g{E). (49) 


Plotting Eq. (55) reveals a Lorentzian-like shape cen¬ 
tered at J = 0. This can also be seen in the ^-dependence 
of the numerical integration of Eq. (54) 


1 (g^ -|- O’) 

En e + ^)^ + ’ 

^ *<!c 


(56) 


As this function has poles when E = —uj, we will now 
integrate around —ui 

/ — LJ+T] 

dEf{E + u;)g{E) (50) 

-LO—r) 

After the change of variable u = E + uj, we obtain 

[ duf{u)g{u - uj) (51) 

J-T) 

The density of states can again be Taylor expanded 
around — w, giving 

g{u-uj) = g{-uj)+ug'{-uj)+ 0{u^). (52) 


The same steps that led us from Eq. (46) to Eq. (48) now 
give us 

2g'{-uj)ri - dg'{-uj)7r, (53) 


which is again linear in d. 

In conclusion, the behavior of the ZPR of non-polar 
materials for non-extremum k-points in the BZ is linear 
when 5 —>■ 0. 


where Sq^ is a sphere of radius gc. This dependence is 
shown at the bottom of Figure for a cutoff radius gc = 
0.5 and uj = 0.01. 

In practice, we also observe that non-adiabatic ZPR 
for VBM/CBM of non-polar materials can be accurately 
fitted by a Lorentzian function. We thus use this type 
of functional dependence of the ZPR in <5 to extrapolate 
the results at i5 = 0. 


2. Polar materials in the non-adiabatic approximation at 
the VBM/CBM 

For polar materials, we get a supplementary ^ factor 
in the integrand of Eq. (54), which thus becomes 

1 


5ft 


1 


d 


q"^ q"^ + U! + id 

= 47r5ft 


dq- 


1 


(57) 


q"^ + UJ + id 
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which becomes 


2tt^ 



dq 


1 

q'^ + UJ + i6 


(58) 


This integral can be performed by closing the contour 
of integration using a half-circle of infinite radius in the 
upper complex plane (which does not contribute to the 
integral) and then using the residue theorem. We obtain 


27r3? 


dq- 


iS 


= 27r^ 




(59) 


which gives 
27r2- 


V^(l + (^) 


S_\2\ 4 



/ S s\ 

1 

( tan 

(-) 

V2 



(60) 


Plotting Eq. (60) again reveals a Lorentzian-like shape 
centered at 5 = 0. Accordingly, fitting the results using 
a Lorentzian is also found to be a good approximation in 
practice for polar materials. 


D. Decreasing 5 does not lead to convergence 


The adiabatic ZPR for the VBM/CBM of a polar ma¬ 
terial has already been shown to diverge with increasing 
q-point sampling (see subsection IV C| ). We now examine 
the d-dependence of the ZPR, wl ' ' 


5ft 


d^q 


1 


rich has the form 
1 


q'^iq'^ + i6) 


= 47r5ft^ 

rQc 

/ dq- 
'0 ' 

47r j 

dq 

Vs 

3 Vs 


q^ -b iS 


52 


(61) 


When qclVS is rather large (i.e. > 10), which happens 
for small 6, the integral of Eq. (61) converges logarithmi¬ 


cally to its value at infinity and we obtain 


We therefore see that Eq. (611 will numerically diverge 


as for small values of S. 

VO 


An example of this kind 
of divergence is given in the middle right of Eigure 9 of 
supplemental material^^ for the VBM of /3-AIN (a polar 
material) using the adiabatic equation. 


VI. RESULTS ON DIFFERENT SEMICONDUCTORS 

We will now examine five different semiconductors, two 
of which are non polar materials (diamond and silicon), 
and three of which are polar materials (a-AIN, /3-AlN and 
BN). We will be able to provide fully converged results, 
independent of any arbitrary parameter, like an ad hoc 
broadening parameter. Of course, as outlined previously, 
in the case of polar materials, only the non-adiabatic the¬ 
ory can provide such results, as the standard adiabatic 
AHC theory breaks down for these. 


A. Non polar materials 
1. Diamond 


Diamond is a metastable allotrope of carbon where 
the C atoms are arranged into two interpenetrating face- 
centered cubic lattices shifted along the body diagonal 
by I**' of its length. The space group associated with 
this spatial arrangement is Ed3m (cubic, 227). Diamond 
has the highest hardness and thermal conductivity of any 
bulk materiaP^. It is therefore used as cutting and pol¬ 
ishing tool in the industry. Even though the stable phase 
of bulk carbon under normal condition is graphite, we 
will focus on the diamond phase. 

The pseudopotential was generated using the fhi98PP 
cod^^with a 1.5 atomic unit cut-off radius for pseudiza- 
tion. The valence electrons of carbon, treated explicitly 
in the ab-initio calculations, are the 2s^2p^3d° orbitals. 

Careful convergence studies (error below 0.5 mHa per 
atom on the total energy) led to the use of a 6x6x6 F- 
centered Monkhorst-Pack k-point samplinj^ of the BZ 
and an energy cut-off of 30 Hartree for the truncation 
of the plane wave basis set. The Perdew and Zunger 
parametrization of LDAS^I was used. The relaxed lattice 
parameter is calculated to be 6.652 Bohr, 1.3% below 
the experimental value o f 6.7 40 Bohr, measured at room 
temperatur^^ (see Table IV for more information on the 
structural properties). 

The electronic bandstructure was computed at the 
DPT level and gave a direct bandgap at F of 5.67 eV and 
an indirect F — 0.727X bandgap of 4.25 eV, intrinsically 
below the experimental bandgap of 5.48 eV at 0 KP( see 
Table 0. 

For the calculations of the ZPR, we used 10 bands to 


describe the active space in Eqs. (15) and (16). 


The convergences with respect to q-point integration 
for the band edges and the direct band gap of diamond 
are shown on Figure where the densest used grid is 
a 125x125x125 q-grid (43680 q-points in the irreducible 
Brillouin-Zone (IBZ)). The P^g state is not the bottom 
of the conduction band, and therefore there are other 
states in the BZ with close energy. This leads to nu¬ 
merical instabilities as the denominator of the adiabatic 


Eq. (15) can diverge for small iS. For large enough imag¬ 


inary component, the ZPR converges to approximatively 
-270 meV. The top of the valence band P^s converges 
linearly and can be extrapolated to infinitely dense q- 
grid. To obtain a definite value for the ZPR, we have to 
converge the ZPR for vanishing 5. The convergence can 
be found on Figure (5b I and shows that large q-point 


grid are required to enter the expected linear regime 
(see section VB for more information). The extrapo¬ 
lated ZPR is -277.61 meV. The VBM can be smoothly 
extrapolated using a square-root fit to zero value of 5 to 
give 160.96 meV. The adiabatic direct bandgap ZPR of 
diamond is computed to be -438.6 meV. For the non- 
adiabatic direct bandgap of diamond, the convergence 
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lattice parameters [Bohr] 



Space group 

Ecut [Ha] 

k-grid 

this work (LDA) 

other DFT (LDA) 

other DFT (GGA) 

experiment (300K) 

a-AlN 

PGsmc [186] 

35 

6x6x6 

5.783/9.255 

5.820/9.33^ 

5.913/9.48P' 

5.881/9.41^51 

5.880/9.40^ 

5.877/9.41iI511 

P-AIN 

F43m [216] 

35 

6x6x6 

8.130 

8.20d^ 

8.16#^ 

8.31'P' 

8.30d^ 

8.25d^ 

8.27’^ 

c-BN 

F43m [216] 

35 

8x8x8 

6.746 

6.75#^ 

6.75#^ 

6.81#^ 

6.83#^ 

6.85^^ 

6.83ll^ 

6.83^551 

C-d 

Fd3m [227] 

30 

6x6x6 

6.652 

6.65i^ 

6.75#^ 

6.74cP 

Si 

Fd3m [227] 

20 

6x6x6 

10.170 

10.22^1 

10.33^^ 

10.2d^ 


TABLE IV. Convergence parameters for the different studied compounds. The space groups are given in Hermann-Mauguin 
notation with the number in bracket being the crystallographic index number in international Tables and the homogeneous 
k-points sampling are F-centered. All the pseudopotentials in this work use the LDA exchange-correlation functional. 




direct 

gap [eV] 



indirect 

gap [eV] 


this work 
LDA 

other DFT 

exp. 

this work 
LDA 

other DFT 

exp. 

LDA 

GGA 

LDA 

GGA 

a-AlN 

4.691 

4/^33 

4.05#51 

6.28=®^ 

_ 

_ 

_ 

_ 



4.P51 


6.28^ 

_ 

_ 

_ 

_ 



4 4 Jim 


6.^221 

_ 

_ 

_ 

_ 



45 ^ 



_ 

_ 

_ 

_ 



4 7^311 



_ 

_ 

_ 

_ 

/3-AlN 

4.677 

4 JIT] 

3.99^251 

_ 

3.308 


3.30d^ 

_ 



4jli 


_ 


sjni 


_ 



4.3^ 


_ 


3^46] 


_ 



47 ^ 


_ 




_ 

c-BN 

8.890 

s.fpsi 


14.^ 

4.446 

4.4231 

4.45d22l 

6.#11 



8.is^ 




5.1#3 


gM 



8.^ 







G-d 

5.670 


5.57P^ 

7 34 IIII 

4.250 


4 1134111 

5.48® 








4.12® 


Si 

2.567 

2.5^ 

2.557=**211 

3 378416] 

0.463 

OA^ 

0.612® 



TABLE V. The direct and indirect (when relevant) DFT electronic bandgaps are compared with other references (theoretical 
or experimental). The star * sign denotes low temperature experiment (below lOK) and no star means room temperature. 


can be found on Figure]^ and shows that the Ffg state 
converges similarly but the VBM has a rapid convergence 
in q-point integration and a Lorentzian behavior for the 
6 extrapolation. The fitted Lorentzian have three fitting 
parameters: a multiplicative constant A, the full width 
at half maximum (FWHM) and an additive constant B 

jy 

- +B, (62) 

(5)2 + cc2 ’ ^ ^ 

where here A = 10.11, F = 0.55 and B = 121.90. The 
extrapolated ZPR is of -283.23 meV and 133.57 meV 
for the Fjj and F 25 states, respectively. This leads to 
a reduction of -415.8 meV of the direct bandgap due to 
electron-phonon interaction at 0 K. 

The convergences of the CBM for diamond are given in 
Figures 1 and 2 of the supplemental material^^ using the 
adiabatic and non-adiabatic equations, respectively. For 
the adiabatic case, the fact that the q-convergence is not 


smooth for relatively small iS is due to the finite k -|- q 
sampling. Indeed, when we compute the renormaliza¬ 
tion at one of the 6 symmetry equivalent CBM k-points, 
the k -t- q sampling is such that the other five equivalent 
k-points are not sampled exactly (not the numerically 
accurate minimum). The extrapolated ZPR of the CBM 
state is -219.24 meV using a square root fit for the adi¬ 
abatic equation and -196.22 meV using a Lorentzian fit 
for the non-adiabatic equation (see Table VI for more 
information). 

The temperature dependence of the direct and indirect 
bandgaps is reported on Figure for a 75x75x75 q-grid 
and shows that the slope at high temperature for the non- 
adiabatic renormalization with a Lorentzian extrapola¬ 
tion to vanishing imaginary parameter 6 is -0.504 meV/K 
for the direct bandgap and -0.435 meV/K for the indi¬ 
rect one. The phonon-induced broadening 


2t 


{adiabatic,RI A) 


of Eq. (18) is calculated for the 75x75x75 q-grid to be 
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(a) Adiabatic ZPR of direct-gap diamond. 


(b) iS extrapolation of the ZPR. 


FIG. 5. Convergence study for the adiabatic (a) q-point grid density and (b) i5 parameter for the direct bandgap ZPR of 
diamond. The bottom Figures are the difference of the two Figures above them. The adiabatic ZPR of the direct bandgap of 
diamond is -438.6 meV. 


180 meV and 63 meV for the direct and indirect bandgap 
of diamond at 0 K, respectively. 

Our result underestimates the experimental ZPR of the 
diamond indirect bandgap of -364 meV^by 9.4%. Since 
our calculation neglects several effects like anharmonicity, 
non-rigid-ion terms or many-body GW corrections, we 
are rather close to the experimental value. The measured 
linear slope at high temperature for the indirect bandgap 
of diamond is -0.54 meV /KP^. The measured linear slope 
at high temperature for the direct bandgap is —0.60 or 
—0.69 meV/K!^, depending on the analysis^. Our theo¬ 
retical values for the indirect and direct bandgaps under¬ 
estimate the experimental ones by 19% and 16%, respec¬ 
tively. We hypothesize that this underestimation of the 
linear slope at high temperature for the direct bandgap 
of diamond is due to the underestimation of the ZPR 
within DPT. Indeed, as discussed in Ref. na the cor¬ 
rection brought by GW to the ZPR is quite substantial 
for the direct bandgap (-209 meV). Since the ZPR is di¬ 
rectly linked with the slope at high temperature, it is not 
surprising that we witness such an underestimation with 
respect to the experimental results. We did not compute 
such GW correction for ZPR of the indirect bandgap of 
diamond but expect from the results of Figure to have 
a smaller correction. 

It is worthwhile to note that the complete lack of ex¬ 
perimental data for low temperature (T < 100 K) and 
the relatively large error bars (up to ±10 meV) between 


200 K and 350 K generate an uncertainty of several meV 
on the experimental ZPB!^. This calls for new, reliable 
and wide range temperature measurement of the evolu¬ 
tion of the bandgap with temperature in diamond. We 
hope that our theoretical study will stimulate such ex¬ 
perimental interest. 

Finally, the non-adiabatically renormalized electronic 
bandstructure of diamond at OK along the L — F — X 
high symmetry line is shown on Figure for a 75x75x75 
q-point grid for 6 extrapolated to zero linearly and with 
a Lorentzian for the VBM and CBM. 


2. Silicon 

Silicon is a tetravalent metalloid widely used in in¬ 
tegrated circuits and semiconductor electronics. It has 
a diamond cubic crystal structure with a Fd3m (cubic, 
227) space group. Pure silicon is usually not found in 
nature but manufactured to get monocrystalline silicon 
for use in computer microchips. Although the manu¬ 
facturing process is rather expensive, this material is so 
important that about 50.000 metric tons are produced 
per year worldwid^^. 

The pseudopotential used for silicon was generated us¬ 
ing the fhi98PP codd^ with a 1.0247 atomic unit cut¬ 
off radius for pseudization. The pseudopotential is a 
Troullier-Martins with the Perdew/Wang^parametriza- 
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(a) Non-adiabatic ZPR of direct-gap diamond. 
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(b) iS extrapolation of the ZPR. 


FIG. 6. Convergence study for the non-adiabatic (a) q-point grid density and (b) iS parameter for the direct bandgap ZPR of 
diamond. The bottom Figures are the difference of the two Figures above them. The non-adiabatic ZPR of the direct bandgap 
of diamond is -415.8 meV.) 



FIG. 7. Temperature dependence of the diamond gaps using 
the non-adiabatic temperature dependence on a 75x75x75 q- 
grid with Lorentzian extrapolation to vanishing imaginary pa¬ 
rameter S. The slopes at high temperature are -0.504 meV/K 
for the direct gap of diamond and -0.435 meV/K for the indi¬ 
rect one. The experimental points from Clark et al.^and Lo- 
gothetidis et (using first or second-derivative line-shape 
analysis) are shifted so that the lowest temperature point 
matches the theoretical line. 


tion of LDA. The valence electrons of silicon, treated ex¬ 
plicitly in the ab-initio calculations, are the 3s^3p^ or¬ 
bitals. 

Careful convergence checks (error below 0.5 mHa per 
atom on the total energy) lead to the use of a 6x6x6 F- 


centered Monkhorst-Pack k-point samplinj^ of the BZ 
and an energy cut-off of 20 Hartree for the truncation of 
the plane wave basis set. The relaxed lattice parameter is 
calculated to be 10.170 Bohr, 0.9% below the experimen- 
tal value, measured at room temperatur^^ (see Table 
for more information on the structural properties). 


IV 


The electronic bandstructure was computed at the 
DFT level and gave a direct gap at F of 2.567 eV and 
an indirect F — 0.848X bandgap of 0.463 eV intrinsically 
below the experimental bandgaps of 3.378 and 1.17 eV at 
~iokISsI for the direct and indirect bandgap, respectively 
(see Table [v|). Many ab-initio simulations have been per¬ 
formed on silicon and give similar values to ours. For 
example. Ref. 15^ got a direct bandgap of 2.52 eV and an 
indirect one of 0.45 eV, also using the Abinit software. 


For the calculations of the ZPR, we used 10 bands to 
describe the active space in Eqs. (151 and (16). The 
convergence with respect to q-point integration for the 
direct bandgaps of silicon is shown on Figure 3 of the 
supplemental material^^ for the adiabatic equation and 
gives a linearly extrapolated ZPR of -6.23 meV for the 
state and a square-root extrapolation of 40.87 meV 
for the ZPR of the VBM, thus leading to a ZPR of - 
47.1 meV for the direct bandgap. The non-adiabatic di¬ 
rect bandgap ZPR shown on Figure 4 of the supplemental 
material^ of the F^^g state is calculated to be -7.36 meV 
and 34.87 meV for the VBM. The non-adiabatic bandgap 
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OK. figure. 

FIG. 8. Electronic bandstructure (plain black line), renor¬ 
malisation (dashed line) and phonon-induced broadening (en¬ 
velope around the dashed line) at OK using the non-adiabatic 
ZPR integrated on a 75x75x75 q-point grid for <5 extrapolated 
to zero for diamond, where the dots are the actual renormal¬ 
ization calculation. A spline function is used to connect the 
renormalization dots. 


ZPR is therefore slightly smaller than the adiabatic one 
with a value of -42.1 meV. The densest grid computed 
is a 100x100x100 q-grid (22776 q-points in the IBZ). 
The convergences with respect to the indirect bandgap 
of silicon are shown in Figures 5 and 6 of the supple¬ 
mental materialJ^ for the adiabatic and non-adiabatic 
equations respectively. The ZPR of the first one can be 
extrapolated to -23.28 meV for the CBM and 40.87 meV 
for the VBM, leading to -64.3 meV renormalization of 
the bandgap. The second one can be extrapolated to - 
21.43 meV for the CBM and 34.75 meV for the VBM, 
leading to a smaller -56.2 meV renormalization of the 
bandgap. 

Such values can be compared with those mentioned in 
the recently published paper by Patrick and GiustincP^ 
who obtain values of -57 and -22 meV for the indirect 
and direct bandgap renormalization of silicon using a 
4x4x4 supercell within the AHC framework (adiabatic 
equation). Their results matches ours for the adiabatic 
4x4x4 q-grid (1/fVg = 0.25) with -52 and -29 meV for the 
indirect and direct bandgaps, respectively (as can also be 
seen on Figures 3 and 5 of the supplemental material^^. 

The non-adiabatic temperature dependence of direct 
and indirect T — 0.848X bandgaps integrated over a 
75x75x75 q-grid with a Lorentzian extrapolation to van¬ 
ishing imaginary parameter S is reported on Figure 
and gives slopes at high temperature of —0.147 and 
—0.255 meV/K, respectively. 

The experimental zero-point motion renormalization of 
the silicon indirect bandgap is —62 meV obtained from 
mass derivative of the gap or —64 meV obtained from 



FIG. 9. Temperature dependence of the silicon gaps using 
the adiabatic temperature dependence on a 75x75x75 q-grid 
with a Lorentzian extrapolation to vanishing imaginary pa¬ 
rameter S. The slopes at high temperature are -0.147 meV/K 
for the direct gap of silicon and -0.255 meV/K for the indirect 
one. The experimental data from Ref.^ (black circles) ancP^ 
(blue diamond). 

linear extrapolation to 0 KP^. The m easur ed linear slope 
at high temperature is —0.32 meVThose exper¬ 
imental results are larger than our theoretical ones using 
the non-adiabatic extension to the AHC equations, as ex¬ 
pected from DFT calculations. This is linked with the 
fact that we underestimate the ZPR with respect to GW 
calculations. 

Additionally, we present on Figure the phonon- 
induced broadening of the direct and indirect bandgap 
of silicon with temperature. The direct and indirect 
bandgap broadening at 0 K are computed to be 31 meV 
and 23 meV, respectively. The experimental broadening, 
measured with spectroscopic ellipsometry in Ref. [SHI (red 
dots), is attributed to the broadening of the Ei=A 3 — A); 
direct transition with temperature. In Ref. [SHI it is also 
mentioned that the measured values at higher tempera¬ 
ture (black dots) are difficult to attribute to the broad¬ 
ening of one particular transition because the Ei gap is 
nearly degenerate with the Eg = Since the el¬ 

lipsometry measurement is a spectroscopic measurement, 
it can only probe direct transitions. In consequence, we 
should compare the broadening results (both black and 
red dots) with the silicon red line (direct-gap Eg). 

Finally, the non-adiabatically renormalized electronic 
bandstructure of silicon at 0 K along the L — F — X high 
symmetry line is shown on Figure [TT] for a 75x75x75 q- 
point grid with S extrapolated to zero. 


B. Polar materials 

1. Aluminum Nitride 

Aluminum nitride in the wurtzite structure (a-AlN) 
is one of the widest bandgap nitride semiconductor. It 
has a Pbamc (hexagonal, 186) space group. a-AlN is 
used for high-temperature electronics and opto-electronic 
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FIG. 10. Temperature dependence of the silicon broadening 
of the direct and indirect bandgap using the adiabatic tem¬ 
perature dependence on a 75x75x75 q-grid. The experimental 
data (red and black dots) are from Ref. 1661 See additional dis¬ 
cussion in the text. 


Energy [eV] Energy [eV] 



OK. figure. 

FIG. 11. Electronic bandstructure (plain black line), renor¬ 
malization (dashed line) and phonon-induced broadening (en¬ 
velope around the dashed line) at OK using the non-adiabatic 
ZPR integrated on a 75x75x75 q-point grid for S extrapolated 
to zero for silicon, where the dots are the actual renormaliza¬ 
tion calculation. A spline function is used to connect the 
renormalization dots. 


application^^ with a melting temperature of 3273 k!^. 
The zincblende form of aluminum nitride (/I-AIN) has a 
F43m (cubic, 216) space group and has been reported 
to be experimentally metastabl^^. We will study the 
temperature-dependence properties of these two phases 
using the equations mentioned above to compute the 
ZPR. 

Concerning the numerical details of the calculations, 
the aluminium and nitrogen pseudopotentials were gener¬ 
ated using the fhi98PP cod^^with a 1.0247 atomic unit 
cut-off radius for pseudization and a maximum angular 
channel of Z = 2. Both of them are Troullier-Martins 
pseudopotentials with the Perdew/Wan^^ parametriza- 


tion of LDA. The valence electrons of aluminum and ni¬ 
trogen, treated explicitly in the ab-initio calculations, are 
generated for the 3s^3p^ and 2 s^ 2 p^ configuration, re¬ 
spectively. 

Convergence checks (error below 0.5 mHa per atom on 
the total energy) lead to the use of a 6 x 6 x 6 P-centered 
Monkhorst-Pack k-point samplin^^ of the BZ and an 
energy cut-off of 35 Hartree for the truncation of the 
plane wave basis set. 

The relaxed lattice parameters are calculated to be 
a=5.783 and c=9.255 Bohr for a-AIN and a=8.130 Bohr 
for /3-AlN. These values are at maximum 2.3% below the 
experimental one (see Table IV). 

The electronic bandstructures were computed at the 
LDA level. For a-AlN, there is a direct gap at P of 
4.691 eV. This is well above the 4.056 eV bandgap com- 
puted within the Material’s Project using GGAl^ but 
well in the range of other LDA bandgaps (see Table 0. 
These DFT values naturally und eresti mate the experi¬ 
mental bandgap of 6.28 eV at xhe zincblende 

/3-AlN turns out to have an indirect P — X bandgap of 
3.308 eV, almost identical to the Materials Project value 
(see Table 0. The direct gap at P of 4.677 eV is a bit 
above most the the LDA values reported in the Table. 
We have nonetheless to bear in mind that the 4.2 eV 
LDA direct gap calculation performed in Ref. |3HI is done 
at the experimental lattice parameter. 

For the calculations of the ZPR we used 18 bands to 
describe the active space in Eqs. (15) and ([0. 

The q-point integration for the direct bandgaps of a- 
AIN is shown on Figure 7 of the supplemental material^^ 
for the adiabatic equation and diverges for dense q-grid 
as AIN is a polar material. The non-adiabatic direct 
bandgap ZPR shown on Figure 8 of the supplemental 
material^^ converges linearly with the q-point grid and 
the S behavior can be fitted by a Lorentzian function to 0. 
It gives a ZPR of-183.5 meV for the CBM and 194.2 meV 
for the VBM, leading to a ZPR of the direct bandgap of 
a-AlN of -377.7 meV. The densest grid computed is a 
34x34x34 q-grid (2052 q-points in the IBZ). 

For the same reason as a-AIN, the /3-AlN diverges for 
the adiabatic equation and the divergence is shown on 
Figures 9 and 11 of the supplemental materialJ^. The 
non-adiabatic direct bandgap ZPR shown on Figure 10 
of the supplemental material^^ converges linearly with 
the q-point grid and the S behavior can be fitted by a 
linear or Lorentzian function to 0. It gives a ZPR of - 
187.54 meV for the CBM and 226.08 meV for the VBM, 
thus leading to a ZPR of the direct bandgap of /3-AlN of 
-413.62 meV. The non-adiabatic indirect bandgap ZPR 
of /3-AlN is shown on Figure 12 of the supplemental ma¬ 
terials^ and converges linearly with the q-point grid and 
the S behavior can be fitted by a Lorentzian function to 
0. This results into a ZPR of -108.36 meV for the CBM, 
resulting in a -334.4 meV ZPR of the indirect bandgap of 
/3-AlN (see Table VI for more information). The densest 
computed grid is a 100x100x1000 q-grid ( 22776 q-points 
in the IBZ). 
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FIG. 12. Temperature dependence of the a-AlN gaps us¬ 
ing the non-adiabatic temperature dependence on a 34x34x34 
q-grid with Lorentzian extrapolation to vanishing imaginary 
parameter 5. The slopes at high temperature is -0.772 meV/K 
for the direct gap of a-AlN. The experimental data from 
Ref. [68| (black circles) and |69] (blue diamond). 


The temperature dependence of the three gaps is re¬ 
ported on Figures [T^ and 13 The linear slopes at high 
temperature can be extracted to be -0.772, -0.521 and 
-0.763 meV/K for the direct gap of a-AIN, the indirect 
bandgap of /3-AlN and the direct bandgap of /3-AlN, re¬ 
spectively. 

The phonon-induced broadening is calculated for the 
34x34x34 q-grid to be 117 meV for the direct bandgap 
of a-AlN. The broadening of the direct and indirect 
bandgaps of /3-AlN at 0 K integrated on a 75x75x75 q- 
grid are 118 meV and 108 meV, respectively. 



FIG. 13. Temperature dependence of the /3-AlN gaps using 
the non-adiabatic temperature dependence on a 75x75x75 q- 
grid with Lorentzian extrapolation to vanishing imaginary pa¬ 
rameter S. The slopes at high temperature are -0.763 meV/K 
for the direct gap of /3-AlN and -0.521 meV/K for the indirect 
one. No experimental data were found in literature. 

The experimental ZPR for a-AlN has been obtained 
from linear extrapolation to 0 K of the change of the 
direct bandgap with temperature and yield a value of 
-239 meV with a li near slope at high temperature of - 
0.83 meV/KPlSlMl, in relatively good agreement with 
our -0.772 meV/K value. The obvious disagreement 
with our theoretical value for the direct bandgap ZPR 


(-369 meV versus -239 meV) can be attributed to the 
fact that the experimental data set measured by Brun¬ 
ner et are very scarce and on a narrow temperature 
range (4-298 K). As pointed out by Passleil^ for this 
compounds: “this illustrates the great importance of ex¬ 
tending experimental measurements in wide bandgap ma¬ 
terials far beyond room temperature 


Finally, we show in Figures and the non- 
adiabaticly renormalized electronic bandstructure at 0 K 
along the highest symmetry F — M path for a-AIN and 
along the L — F — X path of the j3 phase of AIN. The 
thickness of the lines is associated with the lifetime of the 
electronic state computed with Eq. (18). 


Energy [eV] Energy [eV] 




(a) q-AIN bandstructure at OK. (b) Zoom of the left-hand side 

figure. 


FIG. 14. Electronic bandstructure (plain black line), renor¬ 
malization (dashed line) and phonon-induced broadening (en¬ 
velope around the dashed line) at OK using the non-adiabatic 
ZPR integrated on a 34x34x34 q-point grid for S extrapolated 
to zero for q-AIN, where the dots are the actual renormaliza¬ 
tion calculation. A spline function is used to connect the 
renormalization dots. 


2. Boron Nitride 

Boron nitride (BN) exists in various crystalline forms. 
Its most stable phase under normal condition is an hexag¬ 
onal layered arrangement. It used as a lubricant and an 
additive to cosmetic products. The cubic boron nitride 
(c-BN) has a zincblende structure and is isoelectronic to 
diamond. It is the second hardest material below di¬ 
amond, but its chemical stability is far superior with 
high thermal conductivity and low dielectric constantP^l. 
Boron nitride is not found in nature and is therefore pro¬ 
duced synthetically from boric acid or boron trioxide. We 
will only study the c-BN polymorph here. 

The boron and nitrogen pseudopotentials were also 
generated using the fhi98PP cod^^with a 1.0247 atomic 
unit cut-off radius for pseudization and a maximum angu¬ 
lar channel oi I = 2. Both of them are Troullier-Martins 
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figure. 


FIG. 15. Electronic bandstructure (plain black line), renor¬ 
malization (dashed line) and phonon-induced broadening (en¬ 
velope around the dashed line) at OK using the non-adiabatic 
ZPR integrated on a 75x75x75 q-point grid for <5 extrapolated 
to zero for /1-AlN, where the dots are the actual renormaliza¬ 
tion calculation. A spline function is used to connect the 
renormalization dots. 


pseudopotential with the Perdew/Wan^^ parametriza- 
tion of LDA. The valence electrons of boron and nitro¬ 
gen, treated explicitly in the ab-initio calculations, are 
the 2s^2p^ and 2s^2p^ orbitals, respectively. 

Careful convergence checks (error below 0.5 mHa per 
atom on the total energy) lead to the use of a 8x8x8 F- 
centered Monkhorst-Pack k-point samplin^^ of the BZ 
and an energy cut-off of 35 Hartree for the truncation of 
the plane wave basis set. 

The relaxed lattice parameter is calculated to be 
6.746 Bohr, 1.3% bel ow the experimental value of 
6.833 BohpSSl (see Table IV for more information on the 


structural properties). 

The electronic bandstructure was computed at the 
DFT level and gave a direct gap at F of 8.890 eV and 
an indirect F — X bandgap of 4.446 eV, intrinsically be¬ 
low the experimental bandgap of 6.4 eV at SOOKp^l (see 
Table [V|) . For the calculations of the ZPR, we used 18 
bands to describe the active space in Eqs. (151 and (161. 

The c-BN is also a polar material and therefore di¬ 
verges for the adiabatic equation as shown on Fig- 
16 and 14 of the supplemental material^^. The 


ures 


non-adiabatic direct bandgap ZPR shown on Figure 13 
of the supplemental material^^ converges linearly with 
the q-point grid and the 6 behavior can be fitted by a 
linear or Lorentzian function to 0 and gives a ZPR of - 
301.48 meV for the F)^ state and 200.5 meV for the VBM, 
thus leading to a ZPR of the direct bandgap of c-BN of 
-502.0 meV. The non-adiabatic indirect bandgap ZPR of 
c-BN is shown on Figure 15 of the supplemental materi- 
alP! and converges linearly with the q-point grid and the 


S behavior can be fitted by a Lorentzian function to 0 and 
gives a ZPR of -205.08 meV for the CBM, thus leading 
to a ZPR of the indirect bandgap of c-BN of -405.6 meV 
(see Table |VI| for more information). The densest grid 
computed is a 100x100x1000 q-grid ( 22776 q-points in 
the IBZ). 

The non-adiabatic temperature dependence of direct 
and indirect bandgaps are reported on Figure [l^ for a 
75x75x75 q-grid with extrapolation to zero <5 and give 
slopes at high temperature of -0.639 and -0.521 meV/K, 
respectively. The phonon-induced broadening is calcu¬ 
lated for the 75x75x75 q-grid to be 315 meV and 136 meV 
for the direct and indirect bandgap of c-BN at 0 K, re¬ 
spectively. 

Finally, the non-adiabatic renormalized electronic 
bandstructure of c-BN at OK along the L — F — X high 
symmetry line is shown on Figure [T8| for a 75x75x75 q- 
point grid with S extrapolated to zero. 


VII. CONCLUSIONS 

After a brief reminder of the theory, we present a so¬ 
lution to the divergence problem due to a residual Born 
effective charge stemming from the finite k-point grid 
in numerical ab-initio calculations. We analyze theoret¬ 
ically the q-point convergence for a polar or non-polar 
material; in the adiabatic or non-adiabatic approxima¬ 
tion; for the renormalization of the band extrema or other 
states. We propose an equivalent analysis for the conver¬ 
gence of the imaginary parameter 5 tending to zero. We 
demonstrate that the adiabatic AHC formalism breaks 
down for polar materials and the non-adiabatic AHC 
formalism should therefore be used for these materials. 
We then propose a systematic procedure to converge the 
zero-point motion renormalization (ZPR) and apply it for 
five semiconductors and insulators (diamond, silicon, the 
a and /3 phase of aluminum nitride and boron nitride). 
For these materials, we present the non-adiabatic renor¬ 
malized electronic bandstructure (at the density func¬ 
tional theory level) due to electron-phonon coupling as 
well as the phonon induced lifetime in the adiabatic limit. 
We also show the temperature dependence of their direct 
and indirect bandgaps and compare them with experi¬ 
ment whenever available. The non-adiabatic ZPR at the 
DFT level systematically underestimates the experimen¬ 
tal results (by less than 10%) except in a-AIN where the 
theoretical value is larger than the experimental one. We 
strongly question the validity of the experimental result 
in this case as the experimental ZPR was obtained by 
linear extrapolation to 0 K on a very limited tempera¬ 
ture range (0-300 K) where the linear regime was not yet 
achieved. We therefore believe that the experimental re¬ 
sult for a-AIN from Ref. [SS] is underestimated. On the 
one hand, we hope that our work will emulate experimen¬ 
tal work on wider temperature range. On the other hand, 
this approach might also be used in the future to compute 
more evolved temperature-dependent properties depend- 
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(a) Adiabatic ZPR of c-BN. 


(b) iS extrapolation of the ZPR. 


FIG. 16. Convergence study for the adiabatic (a) q-point grid density and (b) iS parameter for the direct bandgap ZPR of 
c-BN. The bottom Figures are the difference of the two Figures above them. 




ZPR [meV] 

{dGap/dT)T^^ [meV/K] 

Broadening [meV] 

Compounds Gap 

Adiabatic Non-adiabatic Experimental Non-adiabatic 

Experimental 

Adiabatic limit Experimental 

a-AIN r - r 

- 

-377.7 -23^ 

-0.772 

Q n jSMitilltiSi 

117 

d-AiN r - r 

- 

-413.6 

-0.763 


118 

r-x 

- 

-334.4 

-0.521 


108 

c-BN r - r 

- 

-502.0 

-0.639 


315 

r-x 

- 

-405.6 

-0.521 


136 

c r-r 

-438.6 

-415.8 -32(P1,-45CP 

-0.504 

-0.6(P1, -0.6#^ 

180 

r - 0.727X 

-379.3 

-329.8 -36#^ 

-0.435 

-0.54^ 

63 

Si r-r 

-47.1 

-42.1 

-0.147 


31 

r - 0.848X 

-64.3 

-56.2 -6#21^-6#21 

-0.255 


22 ~3#^ 


TABLE VI. ZPR for the different studied compounds as well as the broadening in the static limit at 0 K. The two experimental 
results for the direct bandgap renormalization of diamond from Ref. [58|are extracted using first or second-derivative line-shape 
analysis. The two experimental ZPR for the indirect bandgap of silicon from Ref. are obtained from mass derivative of the 
gap and from linear extrapolation to 0 K. 


ing on the electronic structure, as e.g. optical properties. 
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FIG. 17. Temperature dependence of the c-BN gaps using 
the non-adiabatic temperature dependence on a 75x75x75 q- 
grid with Lorentzian extrapolation to vanishing imaginary pa¬ 
rameter 5. The slopes at high temperature is -0.639 meV/K 
for the direct gap of c-BN and -0.521 meV/K for the indirect 
one. No experimental data were found in literature. 
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figure. 


FIG. 18. Electronic bandstructure (plain black line), renor¬ 
malization (dashed line) and phonon-induced broadening (en¬ 
velope around the dashed line) at OK using the non-adiabatic 
ZPR integrated on a 75x75x75 q-point grid for <5 extrapolated 
to zero for /3-AlN, where the dots are the actual renormaliza¬ 
tion calculation. A spline function is used to connect the 
renormalization dots. 


Appendix A: Derivation of the renormalization factor 

1. Minimization of the variational second-order electronic 
energy 

The variational second-order electronic energy without 
non-linear core correction can be written as (see Eq. (60) 
of Ref. [20]) 


p(2) ^ ^0 

(27r)3 


' BZ 


occ 


(0) 


JO), ,(i) 


k-|-q,k-|-q ^nkl^nk.q/ 


\“nk,qTea;t,k-|-q,kl“nk/ \“nk Tea;t,k,k+q I “nk.q/ 


+ E 


( 0 )| ,( 2 ) 


nk I ^etct.k.k I 


,( 0 ) 


<P{nexc) 


dv? 


n( 0 ) (r) 


|nW(r)|2dr -h 27rno "q ^ 

Gy.^0 I ^ 


+ nW(G = (Al) 


where we follow the notation of Ref. Hi *.e. the super¬ 
scripts (0) and (1) refer to the unperturbed and hrst- 
order perturbation (here in nucleus motion) of the peri¬ 
odic part of the wavefunction, is the spin degeneracy 
factor, jj. is the first-order perturbed potential ex¬ 

ternal to the electronic system that includes the nucleus 
one 


( 1 ) _ ( 1 ) , -( 1 ) 
^ea;i,k-Eq,k “ ^sep,k+q,k + '^ioc,q> 


(A2) 


and £xc is the exchange-correlation energy per electron 


d'^jnexc) 

d'n? 


n(0) 


dvx 


dn 


i(0) 


(A3) 


The second-order change of the nonlocal potential 
■Cggp k.k is given in Eq. (54) of Ref. Hi 
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The first-order change of the local potential of Eq. (23) 
for G = 0 can be written as 




q —>-0 


47r 1 
flo q 




(A4) 


with 


[-Z, + £c, + 0(q^)^. (A5) 


To write Eq. (Al) is a more compact form, we define 


the following vectors 


X,: ^ (G| 


(1) X 
nk,q/ 


Uj A S„Wnk 

A_ /«^l„.(l) 


w 


* = SnWnk 


(A6) 

(A7) 

(AS) 


where Wnk is the weight of the k-point and includes a 
band n dependence to indicate that it is zero for unoccu¬ 
pied states and where the index i stands for the combined 
planewave component G, band index n and wave-vector 
k indices. It will later be useful to note that 


u^x = 


flo 

(27r)3 

fin 


'BZ 


dkJ2 


■s(«ikl«ik,q) 


(A9) 


= f<^HG = 0). 

In addition, we define the following scalars 


A 

a = 




A 

C = 


Stt 1 
flo q^' 
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(27r)3 ^ 


Finally, we introduce the matrix A so that 


ctAx A ^ 
(27r)3 

^nkl^nk.q/ 
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k-rq,k+q 


(f{nexc) 
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n(0) (r) 
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P Q l^q (G)l^ rAm 
2 flo ^ ln-kG|2 ■ 




Using the above definitions (Eqs. 


and (AlO) 


(A13)), Eq. (Al) can be re-written in short-hand nota¬ 


tion 


q = X^Ax -b a(u'^x)(x'^u) 

-b (fox'^u-b xW-f (c.c.))-b c. (A14) 


The physical value of the first-order perturbed peri¬ 
odic part of the wavefunction x is the one that minimizes 
(we will refer to it as xi) 

Jui^q q = ({(5x|(Axi -b (6 -b a(u'fxi))u -b w)} 

-b (c.c.)) = 0. (A15) 

The real part of the quantity between curly bracket {} is 
therefore zero 

3fl{<5x|(Axi + {b + a(u^xi))u + w)} = 0. (A16) 

Eor the preceding relation to hold for any (5x, the quantity 
in parenthesis () must be zero, 

Axi + {b + a(u'^xi))u -b w = 0, 

which leads to 

Xi = —{b + a(u^Xi))A“^u — A~^w. 

We define 

A 


(A17) 
(A18) 
(A19) 
(A20) 
(A21) 
(A22) 
(A23) 

By multiplying the preceding equation by and isolat¬ 
ing u^xib, we obtain 

—6'ufA“^u 

(A24) 


Xi = Xia -bXib, 

where 

A A -1 
Xla = -A W, 

and we are left with 

xit, = -{b + a(u'l'(xiQ -b xi{,)))A”^u. 
By defining 

b' = 6-ba(u'fxiQ), 


Eq. (A21) becomes 

xib = -{b' + a(u'l'xib))A“^u. 


u'I'xib = 


1 -b a(utA~iu)' 


Substituting this result back in Eq. (A23), we obtain Xi;, 

-b' 


Xi6 = 


1 -b a(utA“iu) 


A-"u, 


(A25) 


and, using Eqs. (A19l, (A20), and (A22), we finally have 


A -1 -6-ba(ufA ^w) , 

xi =-A ^w-b^-—^u. 


1 -b a(utA~3u) 


(A26) 


( 2 ) 

value of A^iq q at the minimum xi 


Substituting Eq. (|A17|) into Eq. (|A14|, we obtain the 

(A27) 


ui^q q = 6*U^Xl -b W^Xl -b C. 


Then, substituting Eq. (A26) into Eq. (A27), we finally 
obtain 

, = -wtA-'w + 0+ A-u + (C.O.)) 

1 + aut A-iu 


— Ifopu^A -b a|wf A ^u|- 

1 + aulA'^u 


-. (A28) 
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2. Macroscopic dielectric constant 


From Eq. (B3) of Ref. [20| and Eq. (A9|, we can de¬ 
duce that the second-derivative of the total energy with 


respect to a monochromatic electric field of wavevector 
located inside the first Brillouin zone is 


E 


d( 2 ) 


rin 


(27r)^ 


IBZ 


dky^Sn 

n 

(P{nexc) 


,(i) 


r(0) 




“nk.ql-^k+q.k+q ^rakl^rak.q/ 


(«ik.ql«ik) + (“ikl^ik.q) 


27rr2r 


|nW(G = 0)p 


din? 


l(0)( 


|„(1) ('puz 

|n^^^(r)pdr -|- 27rf2o EH^^=x^Ax -I- x'^u -I- u'^x -I- a(u'^x)(x^u), (A29) 


G^O 


G| 


where we have used the short-hand notation defined be¬ 
fore. In the same spirit as Eq. (A15), we can find the 
value of X that minimizes Eq. (A29) (that we will call 
X 2 ) and deduce 


which gives 


Ax 2 -I- a(u'fx 2 )u -I- u = 0, 


X 2 = (-a(u^X 2 ) - 1)A ^u. 


(A30) 


(A31) 


Multiplying Eq. ( |A30| ) by and isolating X 2 allows us 
to obtain 


X 2 = 


A-^u 

1 -I- o(utA^iu)' 


(A32) 


Substituting Eq. (A30) and then Eq. (A32) into 
Eq. (A29), we obtain the value of at the minimum 


X 2 


pe/( 2 ) _ t„ _ 
^-q,q — U X 2 — - 


1 -I- a(ufA“iu) 


(A33) 


We can also define a total energy where the divergent 
G = 0 Hartree contribution has been removed. The re¬ 
sulting term is analytic 


= x'^Ax -I- x^u u^x, (A34) 

= u'^x -I- x'f(Ax -f- u). (A35) 


The location X 3 of the minimu m of 
tained in a similar way to Eq. (A32) 

Ax 3 


u = 0 ^ X 3 = A ^u. 


(A36) 


value of at the 


Substituting Eq. (|A36| into Eq. (|A34|), we obtain the 

(A37) 


-q.q 


minimum X 3 


^e/^an( 2 ) ^ 


Comparing Eqs. (A33l and (A37), we deduce 


E 


=/( 2 ) 


^ef,an{2) 

*-q.q 


1 -aK 


ef,an(2) ' 

-q.q 


(A38) 


The polarizability x(r, r') is defined as the microscopic 
response to a change of external potential that gives the 
total change of electronic density dn(r) 

5n{v)= ( x(r,r')<5wea;t(r'). (A39) 

'J Sin 

(A40) 


= .(r r') 
(r') 


5v, 


Transforming to reciprocal space and taking the 
(G, G') = ( 0 , 0 ) matrix element yields 


rW(G = 0) = Xq(G = 0,G' = 0), 


(A41) 


where the ( 1 ) superscript refers to the first-order pertur¬ 
bation in the external potential due to the electric field. 
Taking a long wavelength monochromatic electric field as 
the perturbation 


which gives in reciprocal space 


V 


W 

exi ,k+q. 


k(G,G') = <5 g.g' 


Eq. (A33) tells us that 


pe/( 2 ) _ t 

A_q,q - U'X2, 


= ^Xq( 0 , 0 ), 


(A42) 

(A43) 

(A44) 

(A45) 

(A46) 


where the second and third equalities stem from 
Eqs. (A9) and ( |A41[ ), respectively. 

The dielectric function is defined as (see Eq. (23) of 
Ref. [7l] for example) 


HG, G') = <5g.g' + 


47r 


G? 


Xq(G,G'), (A47) 


and the macroscopic dielectric function, which is an aver¬ 
age response to an applied field is (see Eq. (15) of Ref.ITTI 
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for example) 


eM(q) = — 


eq^( 0 , 0 ) 


(A48) 


Therefore, using Eqs. (AlOl, (A46), and (A47), we obtain 


free charge (V • X>(r) = 0 q • X>(q) = 0) and magnetic 
field (V X 5(r) = 0 q X 5(q) = 0) to all ow us to de¬ 
duce the form of the electric field from Eq. (A52) 


_ 47r 

Qot- 

^^0 / Q'y^'ySQS 


(A57) 


eM(q) = 


1 


1 + aE, 


e/( 2 )■ 

-q.q 


Using Eqs. ( A37| ) and (A38) finally yields 


£M(q) = l-aE5q“;(^) 

= 1 + a(u'fA~iu). 


(A49) 


(A50) 

(A51) 


Substituting Eq. (A571 into Eq. (A52) then allows to ob¬ 
tain 


E 


( 2 ) 

-q,q 


= E' 


,an(2) 27r (X /7 


-q,q 


^0 ^2-fS 9'y^'ySQS 

We also introduce a mixed term 

^mix(2) A +A —1 

E __ ' = -w'A u. 


-q.q 


(A58) 


(A59) 


3. Born effective charge 

Eollowing the phenomenological discussion of Born and 
Huang (see p.265 of Ref. [15]), we can extend the total 
energy density Etot (including the vacuum energy) in the 
long-wavelength limit quadratically in ionic displacement 
for the atom k in the direction a and macroscopic 
electric field in the direction a 


Et. 


= b E E S E 

/ X r 

KK 'yd -yd 


(A52) 


k7,(5 
d^E 


Injecting Eqs. (A59), (A55), and (A37) into Eq. ( A28| ), 


( 2 ) 

we can express the the total Alq q as 


/j T^mixKz.) , / \\ 

p( 2 ) ^ pan( 2 ) + (c.C.)) 

-q,q -q,q , f^e/,an(2) 


1 - «^-q,q 


\b\^Ey 


q^q 


1 - 


(A60) 


where the numerator of this equation should be a square 
to establish the connection with the effective charges 
Z*^ ^ of E q. ( A58 ). To make the link, we can compare 
Eqs. (AlO) and (All) to find 




where ilo is the volume, C%^ = ^ is the analytic 

interatomic force constant (IFC),^e„^ = ^ is the di- Ti^g^efore Eq. ^ becomes 
electric function and Z*^ ^ the Born effective 


(A61) 


charge. is associated to the second-order energy (see 

k' jS 


Eq. (All) where the non-analytic terms in q have been 
removed 


4xtAx-f-xW + wtx + 


c. 


(A53) 


^( 2 ) ^^an( 2 )_^| ( 1)|2 

-qq -q,q ^ 


(i.B”f+(«.) 


eM(q) 

I „rmxy^),2 . ‘2 tt |,,,( 

Q|^-q,q I +fUPq 

eM(q) 


(A62) 


The solution X 4 that minimizes Eq. (A53) is 
X 4 = —A^^w, 


(A54) 


which leads to the following value for E“"^'q at the vari¬ 
ational minimum 


= —w'A w + C = G/ 


■an 

KOt . 
KOt 


(A55) 


where we have introduced the macroscopic dielectric 


function using Eqs. (A51| and (A37). The second term 


will be canceled by a contribution from the Ewald ion- 
ion energy to the lowest order in q and can therefore be 
included in the analytic partP^ 


pEw,an(2) _ pan{2) _ _| (1)|2 

-q,q -q,q Un ’ 


(A63) 


The definitions of the electric displacement The non-analytic term in q — 0, i.e. the remainder of 


'D ^ £ + 4:Tt'P and the polarizability 


E^^q q, can be written as 


, _ 1 a(E+^g^) 

' cx. — 


Un 


d£a 


(A56) 


1 ^^-|graW+2Er“('^*|2 


g^eM(q) ^0 


-q,q 


(A64) 


where we have excluded the energy of the electric field 
in vacuum. We can combine this with the absence of 


where we have replaced a and b by their definition 
Eqs. (jAlOj) and (jAllj). 
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For vanishing q we have, at the lowest order (see 
Eq. 1^) 


+c>(q)- 


(A65) 


Moreover, the macroscopic dielectric constant can be 
written as 


1 


£m{<i) = -^^q-y£jsq&- 


(A 66 ) 


7^ 


Eq. (A62) thus becomes 


p(2) p.En,,aui2) 27r ^ 


By identification with Eq. (A58), we deduce 




■mix{2) 

q.q ’ 


(A67) 


(A 68 ) 


where we took the comp lex c onjugate of the quantity 
between the norm of Eq. (A67). 


The total Born effective charge is the sum of the ionic 
charge on the atom k and the electronic charge belonging 
to this atom 




which naturally gives 




797 - 


(A69) 


(A70) 


The last equation leads, in short hand notation, to the 
following relation (see Eq. (|A59[)) 


-u^A ^ E ^Z^a, 




(A71) 


Finally, to the lowest order in q, we deduce from the 
preceding relation, Eq. (|A10) and Eq. (All) 


& —au^A 


= E 

7 

= E 


Airiq^ 

rioq'^ 

Airiq^ 

floq^ 


(Zf^ + (A72) 




(A73) 


4. Derivation of Eq. ([2§ 

The first-order Hartree potential diverges as 1/g be¬ 
cause of a residual electric charge in the first-order den¬ 
sity. The first-order density at G = 0 can be written 
using Eqs. (A9) as 




(A74) 


and (A26l as 


»^q^( 0 ) = 7r( 


rig 


,tA-E 


(&—a(ulA ^w)) 


u^A”^u). (A75) 


1 -I- a(utA“iu) 

Using Eqs. (A51), (A 66 ), ( A71| ), and (A73l, we have 




-yL, 

4772 '7* 


E 47r2 

7 t_A -1. 


o2 1^75 Q'l^'rsqs 


uTA~-"u . (A76) 


Finally using Eq. ( A51|) to replace u^A by —- 

and then Eqs. (lAlOl), (lAfifip, and (|A69|), we deduce 


(0) = - E ^ (^«<5a7 - T 


Z* 


7 ” 9 

The first-order Hartree term 


^75 Q'j^'ysqs 


). (A77) 




( 1 ) 


(0)=4^uE 


2 ’ 


(A78) 


can then be renormalized to account for the slow k-point 
convergence of the Born effective charges by enforcing 
effective charge neutrality within the primitive cell. To 
do so, we introduce the average Born effective charge per 
atom 


-i * n.t 


(A79) 


where Nat is the number of atoms in the primitive cell, 
and subtract it from Z* „ 


V 


ren{l) 

//,q 


( 0 ) 


Uog2 V 


Z 


{Zaa^^ Za^) 

^ yi'yS Q'l^'rsqs 


), (A80) 


which finally gives 


V 


ren{l) 

//,q 


(o) = 


,(i) 




.q( 0 ) 




Q'y ( 1 

\ ^ 






r 


(A 81 ) 


Appendix B: Behavior of the ZPR with the shifted parabola 
model 


In this appendix, we study the behavior of the shifted 
parabola energy model presented in Eq. (35) to mimic 
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the ZPR of polar and non-polar materials at points of 
the BZ that are not VBM nor CBM. 

In spherical coordinates, Eq. (351 can be re-expressed 
as 


e(q) = — 2qqo cos 9 (Bl) 

where we have chosen the shift qo along the z-Cartesian 
axis. This function vanishes when q = 0 or g = 2go cos 9. 
The last root is a sphere centered around qo with radius 
qo- 


1. Integration on the spherical shell of poles 


In this section, the ZPR of a polar material is analyzed 
in the case (5 = 0. The set of poles of the ZPR is located 
on a sphere centered around q = qo. We can introduce 
the new variable q = q — qo and express the ZPR as an 
integral on a sphere of radius qc 


dq d9 dcj)— 


sin( 0 ) 


(q + qo)^ q^ - ql 


= 2t: 


dq 


d9 


sin 0 


9 ^ - do Jo dJ + 2 ggo cos 9 + q^ 


The radial part of Eq. (|B2|) can be integrated to 

l9 + 9o| 


rdq 

do Jo 


d^-dl 


In 


Mg - gol^' 


(B2) 


(B3) 


As we would like to understand the behaviour of the 
poles when integrated, we restrict the integral on a small 
spherical shell around q = 0 with qg — A < q < qg + 
A. Expressing everything in terms of a new variable 
u = q — qo, we deduce 


—27r 


/_A uu + 2qo 
that may be rewritten as 


V M 


/.A ^ 

/ duF(u) — \-G(u) 
J-A u 


ln(|u|) 


(B4) 


(B5) 


with 


As I/m and ln(| M|)/M are odd functions of m, the first 
two terms of Eq. (B 8 ) are zero. The first contributing 
terms arise from F'{0) and G"(0) 


F'(0) f du = 2F\0)A 
J-A 

G'(0) / dMln(|M|) = 2G'(0)(AlnA-A), 
J-A 


(B9) 

(BIO) 


which shows that the integral on the spherical shell be¬ 
haves linearly with the width of the shell, as would any 
regular function do. 

In the case of non-polar materials, the 1/q^ prefactors 
is not present, and this makes the derivation easier as no 
angular dependence is present G(m) = 0. The conclusion 
is nonetheless the same 
(1-f -)2 

F{u) = - = F{0) + F'{0)u + Oiu^). (BIl) 

u + 2qo 


The first non-zero term in the integral of Eq. 
linear in 2F'(0)A. 


(B 8 ) is 


2. Integration of the g = 0 pole 


For non-polar material in the non-adiabatic approxi¬ 
mation, at a point different from the CBM or VBM, the 
function that should radially be integrated is 


3? 


dq^ 


1 


(g^ — qqo COS 9 + UJ + i6) 
The last integral leads to 

sin — qqo cos 9 + uj) 


27r 


dqq 


d9- 


((g2 — qqo cos 9 -I- -I- S^) ’ 

which gives the following radial integral 


(B12) 


(B13) 




/ (g^ + 2qqo -I- \ 

V(g2 - 2qqo+uj)‘^ + 6"^ J ' 


(B14) 


When (5 = 0, this function behaves quadratically when 
q tends to 0 because the lowest order of the Taylor ex¬ 
pansion of the logarithm is linear in q. Therefore, the 
integral on a sphere of radius A of this function is 


F(u) = —27rln(|2(;o 
I -k ^ 

G{u) =2n—^. 

u -I- 2qo 


1 -k ^ 

' 1o 

u -I- 2qo 


(B 6 ) 

(B7) 


The F{u) and G(m) functions are analytic within the 
integration range and can be Taylor expanded. Restrict¬ 
ing the expansion to first-order, one gets 



duF{0)- + G{0)^-J^ 
u u 


+ [ dMF'(0)-f G(0)ln(|M|). 

J-A 


(B 8 ) 


cA 


dq- — In 

4go 


(g^ + 2ggo -I- to)^ 
(g 2 - 2qqo + w)" 


= CA^ + 0{A'^). 


(BI5) 

However, for polar materials, the function to integrate 
is similar to a Lindhard function 


1 


4ggo 


In 


-I- 2qqo + -I- 

(g 2 - 2qqo + uj)^ + 6^ 


(B16) 


This function actually tends to a finite value because of 
the linear behavior of the logarithm in q that cancels the 
denominator. Eq. (BI 6 ) for 5 = 0 gives 


1 


4ggo 


In 


(g^ + 2ggo -f 
(g 2 - 2qqo + a;)^ 


= AA + 0{A^). (BI7) 
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3. Integration on a spherical shell around q — Q 


In this section, we focus only on the ^-behavior of 
polar materials in the non-adiabatic framework, where 
Eq. (B16) has to be integrated. 

The integrand has 3 poles: when q = 0, and at the two 
real roots (if any) of — 2qqo + uj, which we call qi and 
q 2 with qi < q 2 - 

Actually, when q — 0, the integrand does not diverge 
when (5 —>■ 0, as shown in the subsection |B 2| We will 
here focus on g = ^2 as the behavior with respect to S in 
g = gi is similar. 

The integration around g 2 is given by 


/*q2+A ^ 

27 t / dq -In 

Jq2-A 4ggo 


(g^ + 2ggo + ujf + 
(g2 - 2ggo + + 5"^ 


(B18) 


We introduce the change of variable u = q — q 2 , and we 
consider <5 <C A <C g 2 


27r 


du- 


1 


■In 


/ ((gi+g2)(2g2))^ 
/_A““4g2go \ {{q 2 - Q'i)(w))^ + 

This integral can be expressed as 

27T 


(B19) 


4g2go V 


^2A In ^4(gi + g 2 )g|) - 2A ln((J^) 


(B20) 


and evaluated as 
27r 


4g2go 


2A 


In (^4(gi + g 2 )g|) - 2Aln((5^) 


S 2A(g2-gi)/ /A^(g 2 -gi) 


<?2 - <7i 


(i„( 


<72 - <7i 


-4 tan 


(52 

-1 / A(g 2 -gi) 

V 6 


+ 1-2 


(B21) 


As A/6 is large, Eq. (B21) reduces to 


27r 

4g2go V 


^2A In (4(gi + q 2 )qi) - 2A ln(5^) 

_ 5 2 A(g 2 -gi) ^^^ / A(g 2 - gi) \ 

q2- q^ 5 \ 6 J 


0 , TT 

-4- 

q 2 - qi 2 


(B22) 


which behaves as 


Cl - C 26 . 


(B23) 


The same reasoning applies when a; = 0. In this case 
the poles are qi = 0 and q 2 = 2qo. The ZPR behavior of 
a non-extremum point in the Brillouin Zone is thus linear 
in the non-adiabatic and in the static case with respect 
to 6. 
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